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Abstract 

A new method for the stabihzation of the sign problem in the Green Function 

Monte Carlo technique is proposed. The method is devised for real lattice 

Hamiltonians and is based on an iterative " stochastic reconfiguration" scheme 

which introduces some bias but allows a stable simulation with constant sign. 

The systematic reduction of this bias is in principle possible. The method is 

applied to the frustrated Ji — J2 Heisenberg model, and tested against exact 

diagonalization data. Evidence of a finite spin gap for J2/J1 >~ 0.4 is found 

in the thermodynamic limit. 
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As well known the Green Function Monte Carlo method (GFMC) allows to obtain the 
exact ground state properties of a many body hamiltonian with a statistical method. One of 
the most severe restriction is that only positive definite Green function GF can be sampled, 
otherwise the method is facing the well known "sign problem". Approximate techniques like 
the fixed node approximation (FN) have been developed to circumvent the sign problem 
but at the very least they cannot be systematically improved to achieve the exact answer 
within statistical errors. This property has severely limited the applications of GFMC to 
fermions and frustrated boson models. In this letter I propose a new approach to stabilize 
the sign problem, the GFMC with stochastic reconfiguration (GFMCSR), which will be 
shortly described below, revisiting also the basic steps of the standard GFMC on a lattice. 



In order to filter out the ground state of a given lattice hamiltonian H the standard 
power method may be applied iteratively : 



where x represents conventionally the index of a complete basis \x > , -ff^'.a; being the 
corresponding matrix elements of the hamiltonian which in the following are assumed real, 
and A is a positive constant that allows the convergence of ipn to the ground state 'ipoix), for 
large n. In numerical calculations of interesting lattice hamiltonians the dimension of the 
basis grows exponentially with the size and the particle number, though the matrix itself 
is very sparse and all its elements Hx/ .j. , for given x, can be generally computed even for 
large system size. In this case an exact application of (|l|) is impossible unless for few steps. 
A way out is to use a stochastic approach , like GFMC ,which is particularly simple on a 
lattice. 

In order to implement stochastically the iteration (|^) the corresponding lattice GF 
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may be decomposed in the following way: 
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where Px',x is a normalized stochastic matrix, hx > is a normahzation constant and the 
matrix s takes into account the sign of the GF. The typical choice is to take px',x = \Gx',x\/bx 
, bx = J2x' \Crx',x\ and Sx\x = ^S^Gx^x, which is identically one if there is no sign problem. 

In the GFMC method the so called "walker" is defined by a weight w and a configuration 
X.. At a given iteration n the walker is assumed to sample statistically the state ipnix) in 
Eq.®, in the sense that the probability Pn{w,x) to have the walker with weight w (not 
restricted to be positive) in a given configuration x satisfies: / dwPn{w,x)w = ipn{x). Then 
the matrix multiplication (|l]) can be implemented statistically , in the precise sense that 
/ dwPn+i{w,x)w = ipn+i(yX), hj the following three steps: 

1. scale the walker weight by bx'. w' = bxW. 

2. select randomly a new configuration x' according to the stochastic matrix p / . 



3. finally multiply the weight of the walker by the sign factor s / : w' —>■ w's i 



In principle the previous Markov process determines, for large n, the ground state of H even 
with a single walker. In practise it is convenient to use a large number M of walkers, which 
I indicate by (wj, Xj) j = 1, - ■ ■ M, shorthand in the following also by vector notations w, x- 



tially to zero as the Markov iteration MI is repeatedly applied and it is basically impossible 
to reach a reasonably large value of n. 

Recently a remarkable progress in GFMC on a lattice was the extension of the FN to this 
case. The method is based on a definition of an effective GF G{/ ^ which is always positive 
definite but yields a good variational estimate of the energy. For later purposes we define 
this effective GF in a slightly different way, by introducing a parameter 7: which allows to 
sample also the negative elements of the GF : 
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If there is sign problem the average walker sign < s > 
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A - H^^^ - (1 + 7)Vsf (x) if x = x' 

where the diagonal sign-Eip contribution is given by 

H , >0,x'^x 
X , X 

For 7 = the usual formulation 0] is recovered, whereas for 7 > the crossing to the 
negative sign region is allowed so that the exact GF can be written as Gx',x = Sx',xGI, ^ 
where s^'^x is finite and non zero and is determined by the ratio with G and G-^ given 
by Eq.(0) and Eq.(^) respectively. The value of the constant 7 necessary to cross the "nodal 
surface" was chosen to be 1/2 in all forthcoming applications. 

In the basic decomposition the stochastic matrix p^'^x = G^, ^/bx and the normaliza- 
tion coefficient bx = I] 3,' .j. are instead determined only by G-^ . 

By omitting the last step w' — ^ wSx',x in the Markov iteration process MI, the state ipn 
is indeed propagated trough the positive GF G^ . The main property used in the following 
is that at any Markov iteration n we can have a statistic knowledge of both the state 
ijjn{x) obtained with the exact GF and of ipl{x) obtained instead with the approximate but 
positive definite one G^ . To this purpose the j^^ walker is defined by two weights Wj and 
Wj corresponding to the propagation of the walker by G-^ and G respectively. These weights 
act on the same configuration xj. Hereafter the vector w_ represents therefore a shorthand 
notation for the 2M components Wj,Wj for j = 1, . . . M. 

The walker vector w_, x allows to determine statistically the state: 

Ipnix) = / d[w] E Pn{m,x) ^Sx,x,Wj/M (6) 
^ j 

and analogously ipl{x) by replacing the weights Wj with the positive ones Wj in the previous 
equation. In this way the configurations generated by the described Markov process MI, 
if weighted with the constants Wj, are distributed for large n, according to the variational 



state corresponding to . This is a reasonable variational wavefunction (WF), which will 
be the initial approximation to which systematic corrections will be applied, as described 
later on. 

Apart for the previous technical definitions, we can explain in few words the basic idea 
used for the stabilization of the sign problem. The iteration MI converges to the ground 
state, but due to the sign problem, only few iterations can be performed with a reasonable 
statistical accuracy. However, the representation of the state ipn{x) in terms of the walker 
population Xj,Wj is not unique. In fact it is perfectly possible to represent the same state 
ipn{x) either with a walker population with very small average sign or with a one with a very 
large average sign. If such reconfigurations are possible each few kp steps, the average sign 
may be stabilized to a large value during the iteration (P and there will be no difficulty to 
sample the ground state for n oo, with no sign problem. 

1 will show that this reconfiguration is well defined and indeed possible. The set of M 
walkers {w,x) are defined via their probability function Pn{w,x) which in turn defines the 
state ipnix) by Eq.(|^). The task is to change P„ onto a new probability distribution 
corresponding to a steadily high sign for the walker population. This without changing the 
information content, the state ipnix). 

Let us define the new state V'n(^)? the one obtained by averaging over in Eq.(|), 
then the reconfiguration is exact if P^ is such that: 

ipn{x) = ipn{x) for all X (7) 

In general it is difficult or impractical to realize all these conditions (0) as their number 
equals the dimension of the Hilbert space. I consider therefore a set of operators O^, k = 
1, ■ ■ ■ p << M and require only p + 1 stochastic reconfiguration conditions: 

EO.^,.V^;(a:) = EO.^>„(x) (8) 

for = 1, ■ ■ - p, beyond the normalization one J2x i^'{x) = J2x i^n{x) 

The previous equations (H) mean that the so called "mixed averages" of the operators 
O'' coincide before and after the reconfiguration. ^ 
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The main idea of this work is that these p + 1 conditions can be fulfilled exactly (for 
chosen operators) by defining the reconfiguration in the following form: 

r M (Ej\Px,\S > y,yj. ] 



(9) 



/ Px ■ 

where f3 = t-^, — ^ is the average sign after the reconfiguration which is supposed to be much 
higher to stabilize the process. The new configurations are taken randomly among the 
old ones {xj}, according to the table p^^, defined below.. The positive weights Wj represent 
a good starting point for the definition of a reconfiguration with large (3. Though there is 
some arbitrariness in the definition of the coefficients pxj , a simple and convenient choice is: 

k 

— f 

where = J, ^ / are the averages over the positive weights of the mixed estimates 
= 0^',x-i corresponding to the operator O'' and the configuration Xj. 
Then, in order to satisfy the WF conditions (||), by using the definition (|^), it is sufficient 

that the coefficients p^ satisfy the following Markovian conditions: 
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which in turn determine the unknown variables a^, for k = 1, ■ ■ -p, for given w,x- 

For hamiltonian not affected by the sign problem [G^^ = G = and f3 = 1) this 
reconfiguration was already used to control the walker population size without introducing 
any source of systematic error. The present more general reconfiguration (^) can be easily 
and efficiently implemented in a similar way. 

Obviously the reconfiguration conditions (H) are equivalent to the exact ones (|^), when the 
number p of linearly independent operators considered in (^) is equal to the large dimension 
of the Hilbert space. An important applicative issue is whether GFMCSR converges, within 
a reasonable accuracy, even with a small number p of meaningful operators O^. 
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We consider the frustrated Ji — J2 Heisenberg spin 1/2 model on a finite square lattice 
with L sites and periodic boundary conditions (tilted by 45 degrees for the L = 32 size only). 
The model hamiltonian is determined by an antiferromagnetic coupling Ji > between 
nearest neighbor spins and a frustrating coupling J2 > between next neighbor ones. [§-[10 



In all forthcoming examples the stochastic reconfigurations were applied frequently enough 
to maintain the average sign before reconfiguration ~ 0.8, condition that minimize the 
statistical fiuctuations. Moreover in each simulation it is important to work with a fairly 
large number of walkers, since in the M — > 00 limit, the GFMCSR results are practically 
independent of the frequency of reconfigurations, as well as the overall constant energy shift 
A. 

The accuracy of GFMCSR for the ground state is displayed in Tabjl] , and compared 
with other methods. The variational WF (used also for GFMC importance sampling 
contains a Jastrow like factor 

^ R,R' 

to mimic the interaction between the spins Sf^ = ±1/2 at sites R, R', where 77 is a variational 
parameter and the two-spin interaction v can be derived by using the method described in 
pT| , yielding an explicit Fourier transform for v. 



Vg/2 = 1 - 



2 — (t(1 — cos qx cos qy) + cos qx + cos qy 



^ 2 — (j(l — cos qx cos qy) — cos qx — cos qy 
with a = 2J2/ Ji- This potential is not defined for J2/ Ji = 1/2, and in such case I have chosen 
to work with a = 0.8. Restriction to any subspace of total spin projection S^^^ = J2r Sfi 
allows to evaluate the spin gap by performing two simulations for S'f^^ = and S^^^ = 1. 
Henceforth I will use the the same potential v in both subspaces, by optimizing rj for the 
S!ot = energy. 

As shown in the table the accuracy of the variational WF is rather poor, and is con- 
siderably improved by the FN, at least for small J2- This kind of accuracy is however not 
enough to determine the rapid increase of the spin gap as J2/ Ji approaches the value 1/2 of 



the classical transition. Instead, as shown in Fig.(|l|) the GFMCSR allows to achieve a good 
accuracy also on this delicate quantity by considering in the reconfigurations only the en- 
ergy and the spin structure factor = J2r,r' e.'^'^^^~^'^ S^S^, symmetrized over all directions 
and for all non equivalent wavevectors q. Remarkably also mixed averages of correlation 
functions that are not included in such reconfiguration conditions are also significantly 
improved (see table). 

The way GFMCSR reaches the large n limit (at fixed number of operators p) is displayed 
in Fig.(^ where the initial n = distribution was obtained by the FN for 7 = 0. For fixed p 
the algorithm is Markovian and reaches an equilibrium distribution for n — > 00, independent 
of the initial one (see example in Fig.^ where p was changed at the iteration indicated by the 
arrow), this in turn will converge to the ground state distribution for large p. A comparison 
with the standard "release nodes" estimate is also shown in the picture. It is clear that 
there is no hope to obtain meaningful results in this case by the direct sampling of the sign. 
On the contrary this method looks very stable and, though approximate, a convergence to 
a reasonable accuracy is obtained even with a very small number of operators, compared to 
the dimension of the Hilbert space. 

The data shown in the table and in the picture indicate that the accuracy of GFMCSR 
may become rather size independent with a relatively small increase of the operator number 
p. The error to work at finite small p is systematic. Thus there is a considerable cancellation 
of this error for the determination of the spin gap displayed in Fig.(|l]). 

The calculation was therefore extended to the large size system up to L = 100 where exact 
diagonalization is not possible. The spin gap as a function of the system size is displayed 
in Fig.(|^). This figure is consistent with the opening of a finite spin gap for J2/J1 >~ 0.4. 
This gap is certainly not an artifact of the variational WF, which is obviously gapless, as 
also confirmed numerically in the same figure. The present numerical results confirm that 
the transition to a spin liquid state with a finite spin gap but no classical order parameter 
should be close to J2/J1 = 0.4. [|Tg 

This work was supported in part by INFM (PRA HTSC) and CINECA grant. 
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FIGURES 
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FIG. 1. Estimate of the spin gap for several methods: variational (empty triangles), FN 
(empty squares), GFMCSR p=l (empty dots) , GFMCSR (full dots) as in the table for L = 16 
(upper points) and L = 32 (lower ones). The exact results are connected by continuous lines. 
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FIG. 2. Energy per site vs. n for GFMCSR with p = 1 (upper curve to the left of the arrow) 
and p = 9 (remaining curves). The triangles represent the standard method with sign problem, 
i.e. with large error bars already for n > 15. 
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FIG. 3. Size scaling of the spin gap. The dashed hnes are hnear fit of the GFMCSR data 
with p = 9, 14, 20 for L = 36, 64, 100 respectively. Lower curves are the variational estimates and 
continuous lines are guides to the eye. 
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TABLES 



J2/J1 


L 


77 


% VMC 


% FN 


% SRe 


% SR 


0.1 


16 


1.2 


2.84 (2.2) 


0.17 (0.1) 


-0.03 (0.0) 


0.02 (0.0) 


0.2 


16 


1.15 


2.80 (2.5) 


0.41 (0.4) 


0.00 (0.2) 


0.03 (0.0) 


0.3 


16 


1.1 


3.25 (2.5) 


0.87 (0.7) 


0.12 (0.8) 


0.05 (0.1) 


0.4 


16 


0.8 


3.38 (2.4) 


1.76 (3.2) 


0.56 (4.5) 


0.26 (0.2) 


0.5 


16 


0.85 


5.65 (10.9) 


3.84 (8.9) 


2.08 (8.9) 


0.66 (1.1) 


0.1 


32 


1 


1.55 (2.5) 


0.22 (0.3) 


0.05 (0.1) 


0.02 (0.0) 


0.2 


32 


1 


1.78 (2.5) 


0.48 (0.6) 


0.15 (0.6) 


0.05 (0.1) 


0.3 


32 


1 


2.23 (2.1) 


0.85 (0.91) 


0.30 (1.4) 


0.10 (0.0) 


0.4 


32 


0.8 


3.07 (4.0) 


1.61 (3.1) 


0.26 (5.6) 


0.21 (0.1) 


0.5 


32 


0.9 


4.51 (10.0) 


2.92 (7.2) 


1.52 (7.7) 


0.46 (0.9) 


0.1 


36 


1.1 


1.86 (2.8) 


0.21 (0.2) 


0.1 (0.12) 


0.02 (0.1) 


0.2 


36 


1.1 


2.22 (2.8) 


0.47 (0.5) 


0.16 (0.5) 


0.07 (0.1) 


0.3 


36 


1 


2.31 (2.8) 


0.91 (1.4) 


0.35 (2.0) 


0.11 (0.1) 


0.4 


36 


0.8 


3.34 (5.5) 


1.74 (4.5) 


0.51 (6.8) 


0.26 (0.3) 


0.5 


36 


0.9 


5.09 (14.4) 


3.34 (11.1) 


1.83 (11.8) 


0.62 (2.1) 



TABLE I. Percentage error of the energy (square antiferromagnetic order parameter as in 
0) for the various methods: variational (VMC), fixed node (FN) ,p = l GFMCSR (SRe) with the 
energy afone and p = 5, 8, 9 GFMCSR estimate (SR) with the energy and for L = 16, 32, 36. 



The statistical errors are about one place in the last digit. 
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